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Abstract 

In multiscale, multi-physics applications, there is an increasing need for coupling numerical solvers 
that are each applied to a different part of the problem. Here we consider the case of coupling a Lattice 
Boltzmann fluid model and a Finite Difference Navier-Stokes solver. The coupling is implemented so 
that the entire computational domain can be divided in two regions, with the FD solver running on one 
of them and the LB one on the other. 

We show how the various physical quantities of the two approaches should be related to ensure a 
smooth transition at the interface between the regions. We demonstrate the feasibility of the method on 
the Poiseuille flow, where the LB and FD schemes are used on adjacent sub-domains. 

The same idea can be also developed to couple LB models with Finite Volumes, or Finite Elements 
calculations. 

The motivation for developing such a type of coupling is that, depending on the geometry of the flow, 
one technique can be more efficient, less memory consuming, or physically more appropriate than the 
other in some regions {e.g. near the boundaries), whereas the converse is true for other parts of the same 
system. We can also imagine that a given system solved, say by FD, can be augmented in some spatial 
regions with a new physical process that is better treated by a LB model. Our approach allows us to 
only modify the concerned region without altering the rest of the computation. 

1 Introduction 

When it comes to the numerical analysis of fluid flows, one has the choice between many different models. 
On one hand there are solvers based on the discretization of the Navier-Stokes equations, for example by a 
numerical scheme of finite differences, finite elements or finite volumes. On the other hand, a new category 
of solvers have emerged over the past decades that are based on kinetic theories. They describe the fluid 
dynamics at a molecular level and can be seen as discretization schemes for the Boltzmann equation. 

In this paper we consider the possibility of solving separate spatial regions of a simulations with a different 
solver. In particular, we are interested in coupling a finite difference (FD) solver of the two-dimensional 
Navier-Stokes equations with a lattice Boltzmann (LB) method, i.e. a solver for the Boltzmann equation. 
The motivation is that depending on the nature of the region, optimal efficiency may be reached with a 
different solver. The same reasoning applies to the implementation of boundary conditions that are more or 
less naturaly formulated in a given numerical scheme. 

As an example, consider the computation of the drag force experienced by a rigid body placed in a 
uniform stream. During the numerical treatment of this problem, two kinds of boundary conditions need to 
be implemented, one for the boundaries of the body and one for the boundaries of the computational domain. 
In the first case, the emphasis is put on a correct representation of physical properties of the body. This part 
of the domain might therefore preferrably be solved with a LB model, which delivers a microscopic description 
of the physics. Furthermore, the LB model gives a direct access to physical quantities such as the drag force. 
On the other hand, the boundaries of the computational domain are required to simulate a space that 
extends infintely in all directions. A lot of computational space can be saved if those boundaries implement 
an appropriate velocity field (,!,)■ Indeed, the velocity field can be shown to be independent of the geometry 
of the body at a large enough distance. It depends only on the drag force and can be computed analytically. 
The implementation of this analytical solution on the boundaries finds a very natural formulation in the FD 
scheme, which is based on macroscopic variables. 
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In order to couple a LB with a FD model, it is crucial to understand how the LB set of variables is related 
to the FD set and vice versa. Our method follows the same arguments as the ones developed in 2 , where 
the coupling between a LB and a FD solver for the heat equation is presented. The connection between the 
two models is achieved through a first-order expansion of the LB variables around a local equilibrium term. 
One finds that the zeroth-order terms of the expansion are related to the macroscopic quantities (ie.the FD 
variables), whereas the first-order terms depend on gradients of those quantities. 

The paper is organized as follows. Section [5] gives a brief overview of the chosen LB and FD models. 
Conceptual differences between the models are pointed out in view of the formulation of the coupling algo- 
rithm which is develoiped in section|31 This rather technical section contains the first order expansion of the 
LB variables and a guideline to the LB/FD coupling. The algorithm is validated in section0]on a Poiseuille 
flow and is found to be of second order in the velocities. Section [3] draws a conclusion and presents some 
plans for future work. 

2 The numerical models 

The FD model used in our study implements a finite difference scheme on a staggered grid. It is explicit 
in the velocities and implicit in the pressure. The chosen LB model implements the LBGK formulation, in 
which the dynamics are expressed in form of a relaxation towards a local equilibrium term. 

In the following, we will restrict the discussion to an overview and a couple of technical aspects of those 
two models. The interested reader will find more details on the LB model in the references 0, 0] and [H]. 
The FD model is explained in many details in the reference This book offers among others a reference 
to a complete implementation of the FD code in the C language. 

2.1 Spatial discretization 

We consider a two-dimensional, rectangular region 



on which we introduce a grid. This grid is divided into i max cells of equal size in the x-direction and j, 
cells in the y-direction, resulting in grid lines spaced at a distance 



The FD model is based on three quantities that are defined on each cell: the pressure (p), the x- 
component (u) and the y-component (v) of the velocity. They are however placed on a staggered grid. A 
given index of the cell is assigned to the pressure at the cell center, to the x-component of the velocity 
at the right edge and the y-component at the upper edge (cf. Figure The reason for this staggered 
arrangement is that it prevents possible pressure oscillations which could occur had all three variables u, v 
and p be evaluated at the same grid points. 

The LB model uses nine variables k — • • • 8 which are all evaluated at the same location of a cell. 
We fixed our choice on the upper right corner. This is to ensure that the LB and the FD model have a 
compatible interpretation of the location of the domain boundary Stt. Indeed, this boundary is defined on 
a cell edge in the FD model. Considering that most implementations of LB boundary conditions set the 
domain boundary on top of a LB node, this leads us to placing the LB node on the intersection of two cell 
edges. 

The situation is depicted on Figure for a system of extent i max = jmax = 3. It shows also that as a 
result of the staggered arrangement, not all extremal grid points of the FD set of variables come to lie on the 
domain boundary. For this reason, an extra boundary strip of grid cells is introduced, so that the boundary 
conditions may be found by linear interpolations between the nearest grid points on cither side. 



n = [0, 1] x [0, h] G R 2 



Sx = l/i 



max 



and Sy = h/j, 



max 



2.2 The FD model 



The FD model is based on a discretization of the incompressible Navier-Stokes equations 




(1) 
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Figure 1: Choice of indices for FD and LB variables on a chosen grid cell 
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Figure 2: Computational grid representing a domain O of size (i max ■ Sx) x (j max • <5y) with i max = _? max = 3. 
The left hand side depicts the staggered arrangement of the variables over the grid when the domain is 
resolved by a FD scheme. In the case of a LB solver, all variables are located on cell edges, as shown on the 
right hand side. The location of the boundary strip is indicated by a dashed line. 
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and the continuity equation 

divw = (2) 

The computation of the successive iterations (u^,v^\p^) => v^ t+1 \ p( t+1 ^) contains two distinct 

steps: 

1. Resolution of the poisson equation to obtain the new pressure field. This computation utilizes the values 
of the pressure and the velocity at the time t: (u^\ !)W,pC') =>■ (p( t+1 )). In presence of Dirichlet 
boundary conditions, this procedure has a unique solution (except for an integration constant). In 
particular, there is no need for knowing the value of the pressure on the boundary. 

2. Computation of the new velocity field according to a finite difference scheme. It uses the pressure field 
atstept + 1: (««, ( u (*+i), „(*+!)) . 



2.3 The LB model 

The LB model can be interpreted as a discretization of the Boltzmann transport equation on the chosen 
lattice. The possible velocities for the pseudo-particles are the vectors v k . They are chosen so as to match 
the lattice directions: if r is a lattice site, r + v k 8t is also a lattice site. In the present case, we use a so- 
called D2Q9 lattice with nine possible velocities: a zero velocity to describe the population of rest particles, 
four velocities for the horizontal and vertical directions, and four velocities for the diagonal directions. We 
restrict our considerations to a lattice with equal spacing in the x and y directions, i.e. Sx = Sy =: Sr. The 
LB model can also be implemented on different kins of lattices, but it needs to satisfy a couple of isotropy 
requirements. For example, the second- and third-order tensors of the D2Q9 lattice are 

2 

^m k v ka v k i3 = y<W and (3) 
fc 

^2m k v ka v kf 3V k7 = 0, (4) 
k 

where 5 a p is the Kronecker symbol, and b is a constant that is characteristic for the lattice. 

The system is described by 9 corresponding density distribution functions t), k = . . . 8, representing 
the distribution of particles entering site r at time t and moving in direction v k . Physical quantities such as 
the local particle density and velocity are defined from moments of these distributions: 

8 8 

P = ^2 m k.fk and u = ^2m k f k v k . (5) 

fc=0 fe=0 

The m k are the lattice weights: those constant values compensate for the different lengths of diagonal and 
non-diagonal directions. 

Although this fluid model is compressible (the density is space and time dependent), it can be shown 
to solve the incompressible Navier-Stokes equations in the small Mach number regime. The fluid pressure 
is related to the fluid density by the ideal gas state equation p = c? s p, where (? s is the speed of sound. 
Therefore, in the LB model, there is no need to solve the Poisson equation for the pressure. It is all built-in 
in the equation of motion for the f k . However, there exists not always a straightforward way of treating the 
pressure on the boundaries, and an appropriate boundary condition has to be found. 

In the BGK approximation, the particle collision is obtained by a relaxation with coefficient lo to a 
truncated Maxwell-Boltzmann equilibrium distribution function f^ cq \ which depends only on the local fluid 
density and velocity: 

f k (f + Vi8t, t + St)- f k (r, t) = -to (/ fe (r, t) - .^ oq) (r, i)) + ■ v u (6) 

where 

f (eq)/- ,s I , b 1/6 



fk (^>*) = a P y 1 + ~ v iaU a + - y~V la U a J - —U J . (7) 

In these formulae and in the further developments, a repeated greek index inside a multiplicative term implies 
a sum on this index. The term a, like b, is a lattice constant. 
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Figure 3: Subdivision of the computational domain O into a FD subdomain and a LB subdomain (f^)- 
One lattice cell on the interface between 0\ and ^2, at the position i = i; n t, is computed by both methods. 
The boundary nodes of the subdomains are represented by white symbols. They must be implemented by 
means of a coupling term between the two methods. 
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One particularity of this model is that it satisfies the mass and momentum conservation at the molecular 
level: 

^2m k f k (f+ViSt,t + St)-f k (f,t))=0 and (8) 

k 

^2m k v ka (f k (f ' + ViSt,t + St) - f k (r,t)) = St F a . (9) 

k 

3 The coupling algorithm 

3.1 The FD-LB interface 

In this section, we cut O into two subdomains Oi and O2 such that = f2i |J Q2. We apply in f2i the FD 
method and in the LB method respectively. For facility, we assume that the subdomains are rectangular, 
Qi occupying the left hand side and O2 the right hand side of the domain. This procedure can however be 
extended with few changes to a general boundary. 

The way boundary conditions are implemented in the FD and the LB scheme has already been touched 
upon in section^] In particular, Figure |2 describes the position of the extremal grid points (drawn as white 
circles and squares) to which a boundary value must be furnished. On the coupling interface between f2i 
and 0,2, those boundary values are taken from the boundaries of the opposite domain. As a consequence of 
the staggered arrangement of the LB values with respect to the FD values, there is a need for an overlap 
between 0± and 2 - There are several ways the coupling can be implemented. We chose a method in which 
the overlap extends on roughly one lattice site. The position of this site is indexed by i = i; n t (see Figure |3J). 

As a conclusion, the FD domain requires the knowledge of the values for Uij for i = i; nt and j = 1 • • • j ma xi 
and the values for Ujj for i = ij nt + 1 and j = • ■ • j m ax- Those values are easily obtained from the LB field 
that offers a natural access to the macroscopic variables 

The LB domain on the other hand requires the knowledge of the 9 values fk;i,j at i = ij nt — 1 and 
j = 0" - jmax- Those values are more difficult to get at. Clearly, the LB values, offering a description of 
the fluid at a molecular level, contain more information than the FD values. In the next section, we present 
a first-order expansion of the f^ k \ It will show that they depend both of the values of the macroscopic 
variables and their gradients. Furthermore, the results of the first-order expansion will serve as a dictionary 
to convert from FD values to LB values. 

3.2 First-order expansion of the LB equation 

The first-order expansion of the LB dynamics is based on three approximations: 
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1. The streaming operator in equation © is replaced by a first-order time and space series: 

f k (f+ ViSt, t + 5t)- f k (r, t) « St dtf k (f, t) + St (v ia d a ) f k (r, t) (10) 

2. The values fk are split up into an equilibrium an non-equilibrium part. The spatial and temporal 
derivatives of the non-equilibrium part are neglected: 

A:=/i Cq) +4 nCq) and ftif° q) « 0, d Q /< ncq) « 0. (If) 

3. All second order velocity terms are neglected in the equilibrium distribution: 

/r - a /'+ i, p— ( i2 ) 

These approximations are consistent with the point of view taken by a first order Chapman-Enskog expansion 
(see for example 0]). In the following, we simplify the notation and replace the approximation sign (ss) by 
a straight equality sign. 

We start by expressing the conservation formulae © and© at the first order. From JHJ we obtain 

^2 m k (d t fk + v ka d a fk) =0 =^ — + div(pu) = 0. (13) 

k 

We consider the case of an incompressible fluid, thus 

^=div(p^)=0. (14) 

Momentum conservation gives 

^mk(Std t Vkafk + Srv ka Vk(ifk) = => d t (pu a ) + dpYl^ = 0. (15) 

k 

This relation expresses the same physical content as the Navier-Stokes equations We have introduced 
the stress tensor, defined as Tl a {3 = mkVkaVkpfk- A second order analysis of the stress tensor enables to 
verify the equivalence between and (| 1 5(1 . but this point is out of the scope of the present development. 
With those relations, we are ready to consider the expansion of the LB dynamics Using 1(10(1 . we find 

- c/f cq) + b 4±F ■ v k = St d t fk + Sr(c ka d a )f k ^ Std t f[ cq) + Sr(c la d a )f[ cq) . (16) 
The time-derivative of the equilibrium term is further expanded: 

dtf { k q) = d p f { p ] d tP + d pu J { p ] d t {pu a ) n 5 3 ^ (-d p U Q0 + 5tF a ) . (17) 
From the isotropy relations 1(314(1 . we know that 



a 



Thus, 



d Yi a0 = -v z d a p. (18) 
dJ ( cq) = b_^ { _^ Ua0 + StFa) = _ aVkadap+ b_v^ StFa (19) 



Plugging into lfH)|l gives 

St 



,(neq) Ot ( ( V ka U a \\ 

Ik = — [av ka O a p-Vk a O a [ap + bp ^ jj 



= St h - V -^d a {pu p ) (20) 

This equation is the final lead in the chain to convert from the macroscopic variables to the set of LB 
variables. In order to summarize our findings, we remember that the fk are split up into an equilibrium 
and non-equilibrium part: fk = //[ + fk n ° q \ where /^ cq ' ) = fj^^ip, u a ) is obtained from equation 10, and 
^(neq) _ ^( ne( i) (p } d a up) is approximated with the help of equation ((20(1 . 

For some values of the index k, the relation ((20(1 can be further simplified by imposing the continuity 
equation J2J). Table explicits the value of /^ ncq ^ for every k. 
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Table 1: Dictionary for the conversion from macroscopic variables to the non-equilibrium LB term. 
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3.3 Coupling the FD values 

Coupling the boundaries of the FD field to the LB field is a simple exercise. We have seen in section l2~2l 
that the values of the pressure need not be coupled. The two components of the velocity field are computed 
in a straightforward manner from relations (0 . Because the FD and the LB variables are not defined at the 
same point of a lattice cell, we need to adjust the values by a linear interpolation. This leads to the following 
relations: 

«Sj = i/*(«Erf+«w-i) ( 21 ) 

= V2(C +l!i +Ci) (22) 

3.4 Coupling the LB values 

For the computation of the equilibrium terms /^ eq \ we need to obtain the macroscopic variables u, v andp 
from the FD field. This is done by linear interpolation: 







-i,i+i) 


(23) 








(24) 






-lj + l +Pi in ,-lj 


(25) 



A certain difficulty consists in relating the pressure field from the FD simulation to the density field of the 
LB simulation. Indeed, both fields contain a constant additive term which is a priori unknown. Expressing 
this term through an offset p Q { s of the pressure field, one has the following situation: 

„2 ..(LB) 

^ W =P^-Pot s (26) 
Po 

We chose to fix this constant by averaging the FD pressure and the LB density on the interface between Sli 
and f^2 

p=l/h / dr p and p—l/h / drp (27) 

Jon! n da 2 J d^ii n oa.2 

and claiming that the density average is constant on the interface: p = pQ. This leads to 

p = Po ( V —^- + • (28) 



It is not clear if this way of doing is optimal. Another solution might consist in solving locally the Poiseuille 
equation on the boundary of the LB field. This would enable to compute the values for p from LB variables 
only. 

The non-equilibrium terms /^ ncq ' ) are based on the gradients of the velocities. Figure [3] shows that two 
of those gradients can be approximated by a centered difference of half the mesh width: 

'V' J : :., = (29) 

= WC-C-.) ( 30 ) 

One gradient needs to be calculated as a centered difference of mesh width, based on interpolated values of 
the velocity: 

'V-! 1 : ,, = ^k^+i + - (Ch + Cm)). (3D 
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There are not enough grid points at hand for computing the fourth gradient at the same level of precision. 
We are luckily saved by the continuity equation (0) which delivers the requested value: 

fc«E-U = - VE-W ( 32 ) 

Now that all missing variables have been computed, we take a step back and discuss the overall algorithm 
of a LB iteration step. For the purpose of this discussion, the dynamics © are split into two steps. The first, 
the collision step, handles the computation of the equilibrium distribution and maps the "incoming particle 
stream" onto the "outgoing particle stream" /^ out ' ) . It is followed by a streaming step that transports 
the particles by a value of 5t in time and VkSt in space. The details of an iteration step are as follows: 

1. On bulk nodes, pit) and u(t) are computed from the incoming particle densities fl ln \t). On boundary 
nodes, all the values p(t), u(t) and fj^ n \t) are obtained from the variables of the FD field at time t. 

2. All nodes, bulk and boundaries, perform the collision step: fi° (r, t) = (l — Lo)fl ln \r,t)+u>fj, cq \f,t). 

3. The bulk nodes perform the streaming step: fk n \r + VkSt,t + St) — fj, outS> (r, t) , for all r such that 
r + VkSt lies on a bulk node. 

Alternatively, it is possible to extend the streaming step to boundary nodes for those values of that are 
incoming from the bulk of the LB simulation. They are kept unchanged, unlike the remaining set of fj:* 1 ^ 
and the macroscopic variables p and u that are provided by the FD field. In our simulations, this procedure 
seemed to produce results equivalent to those of the proposed algorithm. 



4 Validation on a Poiseuille flow 

We propose a validation of our coupling algorithm on the simulation of a Poiseuille flow. This is a stationary 
flow in a channel of infinite length with no-slip boundaries. The boundaries extend horizontally at a height 
y = and y = L. The fluid velocity is strictly horizontal and does not depend on the x position: v — 
0; d x u = 0. The analytical solution of the Navier-Stokes equations |QJ for this problem is known and 
predicts a parabolic velocity profile: 

u(y) = ±-C(Ly - y 2 ). (33) 

The constant C can be related to the body force (C = f x ) or to the pressure gradient (C = —d x p), depending 
on how the fluid is driven. 

In our example, the fluid is driven by a body force. The left and the right borders implement periodic 
boundary conditions in order to simulate a channel of infinite length. Special care must be taken on specifying 
this boundary in the FD model. Indeed, during the simulation it tends to build up a pressure gradient that 
must be eliminated by imposing a constant value of the pressure on the left and right boundary. 

The simulation is performed on a grid of size « max = 3 and j max = 50. The physical channel widht is set 
to L = 1 and the body force has the value F x = 0.01. The numeric values it s i m are compared to the analytic 
solution it ana from equation l|33[l by means of the overall error (the indices of the formula refer to the LB 
convention introduced on Figure [T|l: 

e= V —/ . (34) 

VEi=o x *w(-,j) 2 

We have run three simulations that can be compared among each other. The first simulation implements 
a pure LB model with bounce-back boundaries, the second simulation a pure FD model, and the third 
simulation is a FD-LB hybrid. In the first case, the no-slip property of the walls is implemented by a 
boundary condition known under the name of bounce-back (see e.g. 0]). In the third case, the top and 
bottom strips of size j' max — 3 are computed by the FD model and the bulk domain by the LB model (see 
Figure 0} . 

A remarkable result of the simulations is that the FD-only model reaches the analytical solution at the 
machine level of precision (10 -15 ). Although there exist LB boundary conditions which obtain the same 
result ([33), their implementation is less natural and straightforward than the one of the FD model. We 
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Figure 4: The computational grid for the simulation of a Poiseuille flow is partitioned into three subdomains 
f2i, ^2 and Q3. The FD scheme is used on the boundary domains f2j and O3, and the LB scheme on the 
bulk domain O2 . 
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Figure 5: Simulation of a body-force driven Poiseuille flow with (1) a LB model, (2) a FD model and (3) a 
FD-LB hybrid (see Figure 0J. The curves show the time-evolution of the error, compared to the analytical 
solution of the Poiseuille flow. 
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further remark that the LB model has a faster convergence (in terms of iteration steps) than the FD model. 
The stationary velocity field of the LB simulation is dock distinct from the analytical prediction, due to the 
limited precision of the boundary condition that is known to be of first order. The hybrid simulation shows 
the expected convergence speed of the LB model and an error due to the limited precision of the coupling. 
However, the error is two orders of magnitude smaller than the one due to the bounce-back boundaries of 
the LB-only simulation. The results of the simulations are shown on Figure 

The order of precision of the coupling can be estimated by varying the grid resolution (j ma x) while keeping 
the physical quantities(L, F x ) constant. Figure H3 plots the error of the stationary velocity field as a function 
of the grid resolution. It appears clearly that the coupling acts like a second-order boundary condition for 
the velocity field. No conclusion can be taken concerning the coupling of the pressure field, because the 
latter is constant in a Poiseuille flow. 

5 Conclusion 

In this work, a LB scheme for 2-D incompressible fluid flows is spatially coupled to a FD scheme on a 
computational domain partitioned in two regions. We present a way to relate the LB distribution functions 
fk with the classical physical quantities and their derivatives. Two particular FD and LB schemes are 
introduced, and a complete coupling algorithm between the two is proposed. At the interface, the LB and 
FD are connected so as to preserve continuity of the physical quantities. The connection between the fk 



L = 1 

Jmax 50 

y 




9 



Figure 6: Error of a FD-LB hybrid Poiseuille flow simulation as a function of the grid resolution. A log-log 
plot shows the coupling of the velocity field to be of second order. 
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variables and the standard macroscopic physical quantities is obtained through the analysis of a first-order 
truncated series around the local equilibrium. The equilibrium part of fk is related to the macroscopic 
quantities and the non-equilibrium part to the gradients thereof. Our coupling methodology is indeed an 
approximation since we neglect higher-order derivatives in the nonequilibrium distributions. A validation 
was performed by simulating a Poiseuille flow with FD boundary strips and LB bulk and comparing it with 
an analytic solution. The simulation shows that in this case, the coupling of the velocity field is of second 
order in the grid resolution. 

We consider the work on the LB-FD coupling to be interesting by its own means, as it expresses the 
conceptual differences between the approaches of those models. In particular, it might be inspiring in 
formulating new kinds of boundary conditions for either model. In this sense, the FD model can take profit 
of the physical point of view taken in the LB approach, whereas LB boundary condition can be inspired by 
the strict mathematical formulations of the FD boundaries. 

The value of our hybrid model in practice needs still to be shown, but we are confident that it will 
prove itself useful in a large class of scientific and engineering problems. We plan to provide first sample 
applications by implementing drag force experiments with LB obstacles and FD domain boundaries. 
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